Physica A, Vol.216, 199-212 (1995) 

Simulation of rotating drum experiments 
using non— circular particles 



Volkhard Buchholtz@0, Thorsten PoschelB! 
and Hans-Jiirgen TillemansS'i 

a Humboldt-Universitat zu Berlin, Institut fur Physik, Unter den Linden 6, 
D- 10099 Berlin, Germany 

b Arbeitsgruppe Nichtlineare Dynamik, Universitdt Potsdam, Am Neuen 
Palais, D-14415 Potsdam, Germany. 



http: / /summa. physik.hu-berlin.de: 80 /^thorsten/ 



Hochstsleistungsrechenzentrum, Forschungszentrum Julich, D-524-25 Jiilich, 
Germany and 

Universitdt zu Koln, Institut fur Theoretische Physik, Ziilpicher Str. 77, 
D-50997 Koln 



Abstract: We investigate the flow of granular material in a rotating cylinder 
numerically using molecular dynamics in two dimensions. The particles are 
described by a new model which allows to simulate geometrically complicated 
shaped grains. The results of the simulation agree significantly better with 
experiments than the results which are based on circular particles. 



1. Introduction 



The flow of granular material in a rotating cylinder or drum reveals many 
interesting effects and has been studied by physicists experimentally as well 
as theoretically in various articles. Moreover, the knowledge of the dynamics 
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of the material is of high interest in engineering since many important techno- 
logical processes are based on the movement of granular material in rotating 
cylinders, such as mixing and milling (see e.g. [0-0]). 

Nakagawa used magnetic resonance imaging (MRI, NMR) to observe mustard 
seeds moving in a rotating drum [|J . Seeds are suited because their oil contains 
free protons which are necessary for this method. The velocity and the concen- 
tration of the seeds in the drum could be determined very elegantly with this 
method. Woodle and Munro || and Rajchenbach |J found experimentally a 
hysteresis for the transition from stick-slip motion of the particles to continu- 
ous motion when changing the rotation velocity of a rotating drum. There are 
two different critical angular velocities which separate the regimes of stick-slip 
and continuous flow. If the rotation frequency is increased the critical angular 
velocity is higher than in the case of decreasing frequency. Rajchenbach found 
experimentally that the relation between the surface current J, which is pro- 
portional to the angular velocity of the rotating drum, and the surface angle 
of the granular medium obeys a power law |J: 

j~(e-e c ) m . (i) 

For the exponent m they measured a value of about 0.5. Rolf and Rothkegel 
used a very interesting and sophisticated experimental setup to examine the 
collision frequency between particles in a rotating drum Some of the 

particles contained equipment for pressure measuring and microelectronic de- 
vices to store the data or transmitters to transmit the data to an external 
computer. Several regions with different collision rates were found. 

Presently there are two different methods available for the numerical simula- 
tions of these effects. 

(i) A method which is based on a model by Visscher and Bolsterli |7j was 
introduced in 1992 by Jullien et al. to examine a shaken granular medium^]. 
Here the dynamics of the particles is calculated separately for each particle. 
In the case of the shaken box the particles are ordered according to their 
distance from the bottom of the box. Then one calculates the motion of the 
particles beginning with the one which has the lowest distance to the bottom 
while all other particles are fixed. The particle moves until it reaches a local 
minimum where it gets stuck. Then one calculates the next particle while all 
others are at fixed positions and so on. This method is very fast but it has the 
disadvantage that inertia and elasticity are neglected and only particles with a 
restitution coefficient of zero can be simulated. For this reason the model has 
been attacked in ref. ||. Nevertheless under certain conditions this model is 
able to reproduce the experimental behavior correctly. With a model which is 
based on the same method Baumann et al. found an interesting phenomenon 
which is also measured experimentally HIT)|-|r2|j : Two sorts of circular particles, 
originally perfectly mixed, will decompose after only one or two rotations of 
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the drum. The smaller disks gather themselves in the center of the drum. 

(ii) There are several models for the simulation of the dynamics of granular 
material using the molecular dynamics concept (e.g. |T3). Most of the authors 
apply a model introduced by Cundall and Strack |14| and Haff and Werner |T5| 



where it is assumed that the particles are ideal circular disks i with radii 
Ri which interact when they touch each other. If the distance between their 
centers fj is smaller than the sum of their radii they feel a repulsive force. To 
be able to capture the loss of energy and to simulate a restitution coefficient 
this movement is damped. Hence the particles feel a force Fy with the normal 
component FjJ and the shear component F$ 



F, = FS-^h + F*.\ (2) 




with 



Fg = Y • (|f, - f 3 \ -Ri- Rj)* - mf ■ lN • w& (3) 
Fjj = -min{mf ls \v s rel \^\F^\) (4) 

y rel — _ y>. ^ 

m i = , 6 

J rrii + nij 

with the Young modulus Y, the damping coefficients 7^ and 75 and the 
Coulomb friction coefficient /i. Eq. ([|) describes the relative velocity of the 
surfaces of the particles at the point of contact and eq. @ gives the effective 
mass. Eq. ([|) takes the Coulomb friction law into account, saying that two 
particles slide on top of each other if the shear force overcomes fi times the 
normal force. Eq. (|3|) comes from the Hertz law |T^] for the force between two 



rigid bodies which are in contact with each other. This model has been applied 
to simulations of rotating drums by several authors. Ristow investigated the 
particle size segregation in a rotating drum in two dimensions [1^,1^1 , Poschel 



and Buchholtz investigated the irregularities of the particle stream |L9 . 



There are some models to simulate more complex shaped grains which are 
composed of circular disks. In the model by Gallas and Sokolowski [2C] the 



grains consist of two disks connected with each other by a stiff bar. Walton 
and Braun used a three dimensional MD simulation with spheres and with 
more complicated two-dimensional particles consisting of four or eight circular 
disks rigidly connected with each other [f2l]| . With this model they examined 
the transition from stationary to sliding and the transition from sliding to 
raining flow in a rotating drum. Poschel and Buchholtz [^2| , p3l describe grains 



built up of five circular disks where one of them is located in the center of 
the grain and four identical disks are at the corners of a square. Each pair 
of neighboring disks is connected by a damped spring. The latter model was 
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applied to the rotating cylinder |19| and it was shown that the simulation 
results agree much better with the experiment than equivalent simulations 
using circular disks. Especially it was shown that one can reproduce stick-slip 
motion and avalanches which was not possible with circles. The inclination 
and the dependency of the inclination on the angular velocity, however, did 
not agree well with the experiment. 



Another model for non-circular grains has been proposed by Mustoe and 
DePorter |J24| . There the boundary geometry of the particles is defined (in 
local coordinates) by 





rii 


\Vi\ 




+ 









-1 = 0. 



(7) 



By changing rii from 2 to oo the shape of the particles varies continuously 
from elliptical to rectangular. They use a complicated algorithm to detect 
collisions between particles. Hogue and Newland investigate the flow of 
granular material on an inclined chute and through a hopper using a model 
where the boundary of the convex particles is given by a polygon with up to 24 
vertices. To detect whether two particles touch each other one has to calculate 
the intersections between each pair of vertices. During collisions energy is 
dissipated according to Stronge's energy dissipation hypothesis [^6] in normal 
direction, whilst Coulomb's friction law models the energy losses in tangential 
direction. 



Using grains which consist of interconnected circles or of particles described by 
eq. (0) it is not possible to simulate particles with sharply formed corners. For 
some effects it seems to be essential to simulate such particles to reproduce the 
experimental observed effects. This point is discussed in detail in p3|j27|j28[| . 
In this paper we present molecular dynamics simulations where the particles 
are simulated using a more sophisticated model. 



2. The model 



The particles in our simulation consist of four triangles which are connected 
by beams of length L = P ■ to simulate a square particle of side length 
P as shown in fig. [I] (left). The ends of the beam are fixed at the center of 
mass point of the connected triangles in the direction in perpendicular to the 
neighboring sides of the triangles when the grain is in its rest state, i.e. no 
forces are applied. The beams are shadowed in fig. (JJ). A beam in our sense 
is an elastic bar which is subjected to forces in the direction of its axis and 
transverse to its axis, i.e. to normal and shear forces, and to torques acting on 
its ends. During its deformation a beam dissipates energy similar to a linearly 
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damped spring, i.e. proportional to its deformation rate. 



P 




Fig. 1. Draft of a complex particle (left). The four triangles are connected by damped 
leaf springs as shown in the right figure. The beams are fixed at the centre of mass 
points of the triangles Ca and Cb- When the triangles undergo collisions the beams 
are elongated and bent. They disspipate energy proportionally to the deformation 
rate. 



The deformation of a beam between the triangles A and B with centre of mass 
points Ca and Cb is completely described by the angles 9a and 6# and by 
its length L (see fig. [I], right). A deformed beam acts on the triangles with the 
torques 

M A = A -^e A + 2 -^e B , (8) 
M B = 2 -fQA + A -fe B . 0) 



and with the forces in perpendicular to CaCb 

Fl= 6 ~§- (&a + &b) , (10) 

n = ~ (&A + &B). (11) 

where I is the moment of inertia of the beam and E is the elastic constant 

of the beam material. In the case of a beam with rectangular cross section of 
side lengths S and unity one finds 

I = —S 3 . (12) 

12 v ; 



In the direction of CaCb the beam acts with the forces 



Fa = E ■ (L — L ) , (13) 
Fb — —F ■ (L — L ) , (14) 
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(15) 



The torques and forces given by eqs. |5|-|r3| are valid for small deformations of 



the beam, i.e. we assume a linear superposition of all forces and torques [29] . 
For a detailed derivation of the formulae see [j30f 



The deformation of the beams is damped proportionally to the deformation 
rates 



M 



(d) 
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v A - C A C B -v B - C A C B 
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v A - C A C B -v B - C A C B 



(16) 
(17) 
(18) 
(19) 



where 7 is the damping coefficient of the beam material. As mentioned above 
the deformation of the beam is in linear approximation completely determined 
by the angles Q A and 0s and the length L. Hence there is no damping force 
acting in shear direction. 

When two triangles of different grains touch each other (Fig. |2|), i.e. if there is 
an overlap between both they feel a force in perpendicular to the line between 
the intersection points S1S2 according to the Poisson hypothesis HHIl - The 
absolute value of the force is assumed to be proportional to the intersection 
area Q times Young modulus Y 

Fa \ = Y Q . (20) 
We want to point out that during the collision of the triangles the energy is 



conserved. The method is described in more detail in |30| 



Adding up all forces and torques described above one finds the net forces and 
torques acting upon each particle. The dynamics of the grains will then be 
calculated by a molecular dynamics scheme (see e.g. |I3]] ). We have chosen a 
Gear predictor corrector method of 5th order ||32|| . A beam model very similar 
to ours is described in detail in 



33j. 



3. The experimental setup 



Throughout our simulations we used square particles of side length P con- 
sisting of four triangles as shown in fig. [I] where P is equally distributed in 
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Fig. 2. According to the Poisson hypothesis the force between two interacting particles 
is directed in perpendicular to the intersection line S\S2- 

the interval P G [0.1 cm; 0.2 cm}. The Young modulus Y of the material is 
Y — 2 • 10 7 g ■ cm -1 ■ sec~ 2 . Each connection between the triangles consists 
of beams with elastic constant of the beam material E = 10 5 g ■ sec~ 2 and 
damping constant 7 = 5g ■ sec~ l . The moment of inertia of the beam was 
/ = 10~ 4 cm 3 . When simulating relatively stiff beams we need a very small 
time step St for the integration of the equations of motion to guarantee nu- 
merical stability, and hence the simulation becomes inefficient. An alternative 
method is to simulate each beam by two identical beams of lower elastic con- 
stant E and using a larger time step for the integration. We found that the 
latter approach works more efficiently The parameters have been chosen to 
give the best agreement of the numerical simulated grain movement with the 
typical behavior of granular media. We generated an animated sequence of 
snapshots of the simulation and compared the behavior of the granular ma- 
terial with what we would expect from a real material. The cylinder had a 
diameter D = 4 cm. The wall of the cylinder consists of particles as shown 
in Fig. |H to simulate a rough damping surface. The outer triangles are fixed 
at the uniformly rotating ring, the inner ones obey the laws derived above 
(eqs. |S|-|2"0|) when they undergo collisions with freely moving particles. 

With the integration step St — 2 • 10~ 5 sec the simulation behaves numerical 
stable, i.e. the results do not vary when enlarging the step a few per cent. 

Fig. [| shows a snapshot of the simulation with the angular velocity Q = 
0.2 sec -1 . The iV = 500 complex particles consist each of 4 triangles con- 
nected by 8 beams. The gray scale codes for the particle velocity, black means 
high velocity. The computations have been performed on a IBM 370/R6000 
workstation. The total CPU time for the results presented in this paper was 
approximately 2100 hours. 
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Fig. 3. The wall of the cylinder was built up of particles consisting of two triangles 
connected by damped beams to simulate a rough damping wall. 




Fig. 4. Snapshot of the simulation with N = 500 complex particles with 
P £ [0.1 cm; 0.2 cm] in a rotating cylinder with diameter D = 4 cm. The angu- 
lar velocity is Q = 0.2 sec -1 . The wall particles have not been drawn. The snapshot 
has been taken during an avalanche. The gray scale codes for the particle velocity. 

4. Results 

In our simulation we started with the angular velocity Q = 0.1 sec -1 . As we 
will show below the flow of the material behaves stick-slip like for this velocity. 
After a relaxation time t r = 1.5 sec we began to measure the interesting values. 
These values have been averaged over half a rotation before increasing the 
angular velocity by AQ = 0.1 sec~ l . After reaching the velocity Q = 1.3 sec -1 
where the flow moves continuously we decreased the angular velocity in the 
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same manner while measuring and averaging the interesting data. Fig. [| shows 
the color encoded velocity distribution of the particles moving downwards on 
the surface of the material. The left figure shows the distribution for Q = 
0.1 sec -1 , the right for Q = 1.3 sec -1 . Each horizontal line displays the values 
for a certain time. Time is increased from bottom to top. Red color encodes 
for high particle velocity, blue encodes for low velocity. The color coding is the 
same in the left and the right part of the figure. In the left part one can see 
that for low angular velocity the particles at the surface of the material move 
in a very irregular manner. One observes long periods encoded by blue color 
where almost no material moves downwards, sharply separated from short red 
regions where the material moves rapidly due to avalanches. We want to call 
this behavior stick-slip motion. In the equivalent figure for higher angular 
velocities one does not observe avalanches. The behavior is not stick-slip like 
but continuous in time and space. 

In the natural as well as in the numerical experiment one observes avalanches. 
The size and the time intervals between avalanches vary irregularly. Unfor- 
tunately our numerical data are far from sufficient to calculate reliable size 
distributions of the avalanches. 

In fig. ^|we have drawn the averaged velocity (v s ) of the particles moving at the 
surface of the granular material over time (bottom figure). Each half rotation 
the angular velocity of the cylinder (top figure) was changed according to 
Q = Q + AQ while accelerating, and Q = Q — AQ while decelerating. The 
fluctuations become smaller since the flow is much less irregularly for larger 
rotation velocity. For the case of low rotation velocity one finds avalanches due 
to the stick-slip characteristics while for high rotation speed the flow becomes 
more homogeneous. 

In the simulation we recorded the time series of the material flow gliding down 
at the surface of the granular material /_>. When we draw the relative standard 
deviation 



Tflo 



^ { ^-t )2 (21) 

of this time series as a function of the rotation velocity Q (fig. ^) we find 
that for Q > Q cr ^0.6 sec' 1 the relative standard deviation <Jfi ow is less than 
0.1. For smaller Q the standard deviation <7fi ow rises dramatically due to the 
transition from the continuous regime to the stick-slip motion. For each value 
Q we find two values <7fi ow , one when accelerating the cylinder and one when 
decelerating. The dashed line in fig. [7] connects the mean values of each of 
these pairs. 

Our acceleration-deceleration schedule was intended to reproduce the hys- 
teresis in the critical transition point Q cr where the characteristics of the flow 
changes suddenly, which can be observed in the experiment 0. When accel- 
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Fig. 5. The view at the surface of the moving granular material for consecutive times 
for = 0.1 sec~ l (left) and fl = 1.3 sec -1 (right). The colors encode for the particle 
velocity at the surface, red color means high velocity, blue means low velocity. Time 
rises from bottom to top. For low angular velocities one finds sharply distinguishable 
regions in space and time of different velocities which corresponds to the stick-slip 
like motion of the particles. The flow for = 1.3 sec -1 is more homogeneous. The 
material flows from left to right in each column. 



erating the cylinder the transition occurs at the velocity f2S r which is higher 
than the transition point Qf when decelerating the cylinder. From the results 
presented so far we cannot observe a sharp transition between the stick-slip 
flow and the continuous regime as it is observed in the experiment ||. Hence 
we are not able to reproduce the hysteresis described above. 

Our system is too small to provide a smooth surface to measure the inclination 
directly. Therefore we applied an indirect method which is explained in 
detail in an earlier paper [19| . In fig. |^ we have drawn the inclination of the 
surface of the material for low angular velocity Q = 0.1 sec -1 (left) and for 
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Fig. 6. The averaged positive particle velocity (v s ) over time (bottom figure). Each 
half rotation the angular velocity of the cylinder (top figure) was changed. For high 
£1 the relative fluctuations become small compared with the relative fluctuations for 
smaller Q. 
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Fig. 7. The relative standard deviation (eq. \2J\) of the time series of the positive 
material flow in the rotating cylinder. For each SI we measure two values o~fi ow ac- 
cording to our acceleration- deceleration schedule. The dashed line connects the mean 
values of each of these pairs. We observe a transition between the two different flow 
regimes, but we are not able to determine a precise transition point. The transition 
takes place at the rotation velocity f2 ~ 0.6 sec^ 1 . 



11 



higher velocity Q = 1.3 seer 1 (right). The values at the time axis correspond 
to those in figs. [| and || In the first case the angle is smaller and varies 
much more irregularly due to avalanches. The mass of the material remains 
constant throughout the simulation. Hence the averaged flow over time in 
positive horizontal direction /_> which is due to material flow at the surface 
equals the flow in negative direction /<_ given by the continuous material 
transport due to the rotation of the cylinder. The net flow Af = /_> — /<_ 
fluctuates around zero. When the material moves stick-slip like in a clockwise 
rotating cylinder we expect that there are relatively rare events (avalanches) 
when the material moves in positive direction and time intervals, where the 
flow on the surface is almost zero, i.e. Af = — /<_ (see also fig. ||]). For higher 
rotation speed one expects that the material moves homogeneously and there 
are only small fluctuations of Af. The impulses in the upper part of fig. § 
show the cumulative positive material flow between each two consecutive times 
ti and when Af = divided by the negative flow /«_, i.e. 

e /=o) 

f™(^2^)=J- J A f dt - ( 22 ) 



t (A/=0) 



These values give a quantitative measure for the irregularity of the flow. One 
can see that for low Q one observes fewer but much bigger distinguishable 
avalanches than for high angular velocity. For low angular velocity Q one finds 
a strong correlation between avalanches (top figure and sudden decreases 
of the inclination (bottom figure). 

The inclination of the surface of the material is a function of the angular 
velocity. Rajchenbach || determined the law experimentally to be 

6 - C ~ Q 2 , (23) 

where Q c is a constant. For each angular velocity Q we measured in the sim- 
ulation the inclination of the surface 0. Fig. ^ shows the inclination of the 
material surface over the angular velocity Q and the function in eq. fl2"3"|). The 
numerical data points are located at both sides of the experimental curve. Al- 
though there are not enough data points to confirm the experimental results, 
we can show at least that the experimental and theoretical results do not con- 
tradict each other. The measured value for the critical angle is C = 22.3°. 
Typical values for the experimentally measured angle of repose lie between 



C = 20° and C = 30°, Bretz et al. for instance found C « 25° [34 . Values 



for the critical angle found in earlier simulation |TJ| are C m 8° for circular 



particles and C fa 18° for non-circular compound particles ||19| . 
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Fig. 8. Time series of the inclination of the material over the time (lower fig- 
ures) for the angular velocities = 0.1 sec -1 (left) and Q = 1.3 sec -1 (right). The 
values at the time axis correspond to the previous figures. In the first case the angle 
varies more irregularly due to avalanches. The upper figures show the material flow 
(explanation in the text). 
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Fig. 9. The inclination of the material surface over the angular velocity fL 
The dotted line displays the function which has been measured experimentally by 
Rajchenbach [6j. 

5. Conclusion 



We presented a model to simulate the translational and rotational movement 
of complex particles in a rotating drum. The particles consist of triangles 
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connected to each other by beams. Energy dissipation is caused by inner de- 
grees of freedom. Our model is able to predict the angle of the surface of the 
granular material with the horizontal line with higher accuracy than other 
models. The angle given by our model is in the same range as experimental 
values and higher than the angle predicted by models using circular disks or 
particles consisting of circles. We think that this advantage of the model is 
caused by the use of particles with the more realistic polygonal shape which 
allows to take steric effects into account. One drawback of our model is its 
slowness which is due to the complexity of the numerical calculation. As also 
observed in experiments and in other simulations we found two flow regimes, 
a stick-slip motion of the grains for low drum frequencies and a continuous 
flow regime for high frequencies. Our data fit well with the power law between 
the surface flow and the surface angle found experimentally by Rajchenbach. 
Unfortunately, they are not yet sufficient to calculate the precise exponent. 
Further optimization and more simulation runs have to be performed. 
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